Introduction à la modélisation statistique bayésienne

Thierry Phénix & Ladislas Nalborczyk
17 mars 2018

Planning

Cours

Cours 1: Introduction à l'inférence bayésienne

Cours 2: Modèle beta-binomial

Cours 3 : Modèle linéaire à 1 prédicteur continu

Cours 4 : Modèle linéaire multivarié

Cours 5: Markov Chain Monte Carlo

Cours 6: Modèle linéaire généralisé

Cours 7: Comparaison de modèles

Cours 8: Prédicteurs catégoriels et Interactions

Cours 9: Modèles multi-niveaux

Cours 10: Data Hackaton

Plan

Objectifs
Modèle Gaussien sans prédicteur
Modèle Gaussien avec 1 prédicteur
Régression polynomiale

Objectif de la séance

Interpréter ce genre de données

plot of chunk unnamed-chunk-1

Construire ce type de modèle

\[ \begin{align} y_{i} &\sim \mathrm{Normal}(\mu_{i},\sigma) \\ \mu_{i}&= \alpha + \beta x_{i} \\ \alpha &\sim \mathrm{Normal}(60,10) \\ \beta &\sim \mathrm{Normal}(0,10) \\ \sigma &\sim \mathrm{HalfCauchy}(0,1) \end{align} \]

Planning

Cours

Cours 1: Introduction à l'inférence bayésienne

Cours 2: Modèle beta-binomial

Cours 3 : Modèle linéaire à 1 prédicteur continu

Cours 4 : Modèle linéaire multivarié

Cours 5: Markov Chain Monte Carlo

Cours 6: Modèle linéaire généralisé

Cours 7: Comparaison de modèles

Cours 8: Prédicteurs catégoriels et Interactions

Cours 9: Modèles multi-niveaux

Cours 10: Data Hackaton

Plan

Objectifs
Modèle Gaussien sans prédicteur
Modèle Gaussien avec 1 prédicteur
Régression polynomiale

Modèle Gaussien sans prédicteur

Chargement des données

d <- read.csv("Data/Howell1.csv")
str(d)

'data.frame':   544 obs. of  5 variables:
 $ X     : int  1 2 3 4 5 6 7 8 9 10 ...
 $ height: num  152 140 137 157 145 ...
 $ weight: num  47.8 36.5 31.9 53 41.3 ...
 $ age   : num  63 63 65 41 51 35 32 27 19 54 ...
 $ male  : int  1 0 0 1 0 1 0 1 0 1 ...

Modèle Gaussien sans prédicteur

Distribution de la hauteur

d2 <- d %>% filter(age >= 18)
head(d2)

  X  height   weight age male
1 1 151.765 47.82561  63    1
2 2 139.700 36.48581  63    0
3 3 136.525 31.86484  65    0
4 4 156.845 53.04191  41    1
5 5 145.415 41.27687  51    0
6 6 163.830 62.99259  35    1

plot of chunk unnamed-chunk-4

Modèle Gaussien sans prédicteur

Distribution de la hauteur

d2 <- d %>% filter(age >= 18)
head(d2)

  X  height   weight age male
1 1 151.765 47.82561  63    1
2 2 139.700 36.48581  63    0
3 3 136.525 31.86484  65    0
4 4 156.845 53.04191  41    1
5 5 145.415 41.27687  51    0
6 6 163.830 62.99259  35    1

plot of chunk unnamed-chunk-6

Distribution normale

On suppose que les hauteurs \( h_{i} \) sont i.i.d.
En l'absence d'information, on fait l'hypothèse que les hauteurs \( h_{i} \) sont distribuées suivant une loi normale.

set.seed(1234)    
tibble(value = rnorm(1000, 155, 10)) %>%
    ggplot(aes(x = value) ) +  xlim(130, 180) +
    geom_histogram(bins = 20, col = "white") +
    theme_bw(base_size = 30)

plot of chunk unnamed-chunk-7

Modèle Gaussien sans prédicteur

Distribution de la hauteur

La likelihhod caractérisant les hauteurs \( ~h_{i}~ \) peut-être représentée par le modèle suivant:

\[ h_{i} \sim \mathrm{Normal}(\mu,\sigma) \]

Elle admet deux paramètres \( ~\mu~ \) et \( ~\sigma \).
On cherche à calculer LA distribution postérieure qui décrit le mieux la répartition des tailles \( ~\{h_i\} \).

\[ \color{purple}{p(\mu,\sigma~|~\{h_i\})} \]

\( \longrightarrow \) On va donc explorer toutes les combinaisons possibles de \( \mu \) et \( \sigma \) et les classer par leurs probabilités respectives.

Notre but est de décrire la distribution a posteriori, qui sera donc une distribution de distributions.

Modèle Gaussien sans prédicteur

Calcul de la distribution postérieure

Le calcul de la distribution postérieure des hauteurs est le produit de la fonction de vraisemblance et de la distribution a priori:

Pour la distribution a priori, on fait l'hypothèse de l'indépendance des paramètres \( ~\mu~ \) et \( ~\sigma~ \). La distribution peut alors être vu comme le produit de la distribution de \( ~\mu~ \) par la distribution de \( ~\sigma~ \).

\[ \color{steelblue}{p(\mu, \sigma) = p(\mu)~p(\sigma)} \]

La fonction de vraisemblance est définie par la fonction de densité de la distribution normale:

\[ \color{orangered}{f(h_i~|~\mu, \sigma)=\frac{1}{\sigma\sqrt{2\pi}} \exp\bigg(-\frac{1}{2}\big(\frac{h_i-\mu~}{\sigma}\big)^{2}\bigg)} \]

Modèle Gaussien sans prédicteur

Prior sur le paramètre \( ~\mu~ \)

\( \color{steelblue}{\mu \sim \mathrm{Normal}(178,20)} \)

tibble(
    x = seq(from = 100, to = 250, by = .1),
    y = dnorm(x, mean = 178, sd = 20) ) %>%
ggplot(aes(x, y)) +
    geom_line( col = "steelblue", lwd = 2) +
    xlab(expression(mu)) + ylab("density") + 
    theme_bw(base_size = 30)

plot of chunk unnamed-chunk-8

Prior sur le paramètre \( ~\sigma~ \)

\( \color{steelblue}{\sigma \sim \mathrm{Uniform}(0,50)} \)

tibble(
    x = seq(from = -10, to = 60, by = .1),
    y = dunif(x, 0, 50)) %>%
ggplot(aes(x, y)) +
    geom_line( col = "steelblue", lwd = 2) +
    xlab(expression(sigma)) + ylab("density") + 
    theme_bw(base_size = 30)

plot of chunk unnamed-chunk-9

Modèle Gaussien sans prédicteur

Distribution a priori produit des paramètre \( ~\mu~ \) et \( ~\sigma~ \)

Calcul 'force brute' méthode grille

prior_p <- function(x, y) {
    dnorm(x, mean = 178, sd = 20)*dunif(y, 0, 50) 
}

n <- 200
mu <- seq(from = 100, to = 250, length.out = n)
sigma <- seq(from = -10, to = 60, length.out = n)

prior <- outer(mu, sigma, prior_p)
prior <- prior / sum(prior)

persp3Drgl(mu, sigma, prior, smooth = FALSE, lighting = TRUE, new = TRUE)

You must enable Javascript to view this page properly.

Modèle Gaussien sans prédicteur

Échantillonner à partir de la distribution a priori

set.seed(1234) 
n <- 1e4
tibble(
    sample_mu    = rnorm(n, mean = 178, sd = 20),
    sample_sigma = runif(n, min = 0, max = 50)) %>% 
    mutate(x     = rnorm(n, mean = sample_mu, sd = sample_sigma)) %>% 

  ggplot(aes(x = x)) +
  geom_density(colour = "steelblue", size = 0) +
  scale_y_continuous(NULL, breaks = NULL) +
  labs(subtitle = expression(paste("Prior predictive distribution for ", italic(h[i]))),
       x= NULL) +
  theme_bw(base_size = 30)

plot of chunk unnamed-chunk-12

Modèle Gaussien sans prédicteur

Le piège de la Likelihood

Comme la probabilité de \( ~h_i~ \), les valeurs de \( ~\mu~ \) et \( ~\sigma~ \) étant fixées

\[ p(h_i = 162.8648~|~\mu, \sigma) \] \[ = 0~~~ \forall \mu, \sigma \]

Comme une fonction de \( ~\mu~ \) et \( ~\sigma~ \), la valeur de \( ~h_i~ \) étant fixées

\[ f(h_i = 162.8648]~|~\mu = 151.23, \sigma = 23.42) \] \[ = 0.01505675 \]

où la \( 34^{ème} \) mesure de hauteur vaut \( 162.8648 \)

plot of chunk unnamed-chunk-13

Modèle Gaussien sans prédicteur

Distribution postérieure

\[ \begin{align*} \color{purple}{p(\mu, \sigma~|~\{h_i\})} & = \frac{\color{steelblue}{\mathrm{Normal}(\mu|178,20)\mathrm{Uniform}(\sigma|0,50)}~\prod_{i} \color{orangered}{\mathrm{Normal}(h_{i}|\mu, \sigma)}} {\color{green}{\int \int \prod_{i} \mathrm{Normal}(h_{i}|\mu,\sigma)\mathrm{Normal}(\mu|178,20)\mathrm{Uniform}(\sigma|0,50)d\mu d\sigma}} \\ \\ & \propto \color{steelblue}{\mathrm{Normal}(\mu|178,20)\mathrm{Uniform}(\sigma|0,50)}~\prod_{i} \color{orangered}{\mathrm{Normal}(h_{i}|\mu, \sigma)} \end{align*} \]

Il s'agit de la même formule vue lors des cours 1 et 2, mais cette fois en considérant qu'il existe plusieurs observations de taille (\( h_{i} \)), et deux paramètres à estimer \( \mu \) et \( \sigma \).

La difficulté est le calcul de la vraissemblance marginale (en vert).
Il faut intégrer sur deux paramètres: \( \mu \) et \( \sigma \).

Remarque : la probabilité a posteriori est proportionnelle au produit de la vraissemblance et du prior.

Modèle Gaussien sans prédicteur

Distribution postérieure : méthode Grille

Construction de la grille

n <- 200
d_grid <- tibble(
  mu    = seq(from = 150, to = 160, length.out = n),
  sigma = seq(from = 6,   to = 9,   length.out = n)) %>% 
  expand(mu, sigma) # build all the combination -> matrix 2D

Fonction qui calcule la likelihood.

grid_function <- function(mu, sigma){
  dnorm(d2$height, mean = mu, sd = sigma, log = T) %>% 
    sum()
}

On travaille en log pour réduire les erreurs d'approximation

Modèle Gaussien sans prédicteur

Distribution postérieure : méthode Grille

d_grid <- d_grid %>% 
  mutate(log_likelihood = map2(mu, sigma, grid_function)) %>% 
  unnest() %>% 
  mutate(prior_mu       = dnorm(mu, mean = 178, sd  = 20, log = T),
         prior_sigma    = dunif(sigma, min = 0,   max = 50, log = T)) %>% 
  mutate(product        = log_likelihood + prior_mu + prior_sigma) %>% 
  mutate(probability    = exp(product - max(product)))
head(d_grid)

# A tibble: 6 x 7
     mu sigma log_likelihood prior_mu prior_sigma product probability
  <dbl> <dbl>          <dbl>    <dbl>       <dbl>   <dbl>       <dbl>
1   150  6            -1350.    -4.89       -3.91  -1359.    1.91e-57
2   150  6.02         -1349.    -4.89       -3.91  -1357.    5.73e-57
3   150  6.03         -1348.    -4.89       -3.91  -1356.    1.69e-56
4   150  6.05         -1346.    -4.89       -3.91  -1355.    4.95e-56
5   150  6.06         -1345.    -4.89       -3.91  -1354.    1.43e-55
6   150  6.08         -1344.    -4.89       -3.91  -1353.    4.07e-55

Modèle Gaussien sans prédicteur

Distribution postérieure : méthode Grille

d_grid %>% 
  ggplot(aes(x = mu, y = sigma, z = probability)) + 
  geom_contour() +
  labs(x = expression(mu),
       y = expression(sigma)) +
  coord_cartesian(xlim = range(d_grid$mu),
                  ylim = range(d_grid$sigma)) +
  theme_bw(base_size = 30) +
  theme(panel.grid = element_blank())

plot of chunk unnamed-chunk-19

Modèle Gaussien sans prédicteur

Distribution postérieure : méthode Grille

d_grid %>% 
  ggplot(aes(x = mu, y = sigma)) + 
  geom_raster(aes(fill = probability)) +
  scale_fill_viridis_c(option = "A") +
  labs(x = expression(mu),
       y = expression(sigma)) +
  coord_cartesian(xlim = range(d_grid$mu),
                  ylim = range(d_grid$sigma)) +
  theme_bw(base_size = 30) +
  theme(panel.grid = element_blank())

Cette méthode est beaucoup trop coûteuse et n'est pas utilisée en pratique

plot of chunk unnamed-chunk-21

Modèle Gaussien sans prédicteur

Échantillonnage à partir de la distribution postérieure

Pour échantillonner avec remise sur une grille, on utilise la fonction dplyr::sample_n()

set.seed(434)
d_grid %>% 
  sample_n(size = 1e4, replace = T, weight = probability)  %>% 

  ggplot(aes(x = mu, y = sigma)) + 
  geom_point(size = 2., alpha = 1/15) +
  scale_fill_viridis_c() +
  labs(x = expression(mu[samples]),
       y = expression(sigma[samples])) +
  theme_bw(base_size = 30) +
  theme(panel.grid = element_blank())

plot of chunk unnamed-chunk-23

Modèle Gaussien sans prédicteur

Calcul des distributions marginales \( ~\mu \) et \( ~\sigma \)

set.seed(434)
d_grid %>% 
  sample_n(size = 1e4, replace = T, weight = probability)  %>% 

  select(mu, sigma) %>% 
  gather() %>% 

  ggplot(aes(x = value)) + 
  geom_density(fill = "steelblue", size = 0) +
  scale_y_continuous(NULL, breaks = NULL) +
  xlab(NULL) +
  theme_bw(base_size = 30) +
  theme(panel.grid = element_blank()) +
  facet_wrap(~key, scales = "free")

Les distributions marginales pour \( ~\mu \) et \( ~\sigma \).

plot of chunk unnamed-chunk-25

Modèle Gaussien sans prédicteur

Calcul des distributions marginales \( ~\mu \) et \( ~\sigma \)

library(tidybayes)

set.seed(434)
d_grid %>% 
  sample_n(size = 1e4, replace = T, weight = probability)  %>% 

  select(mu, sigma) %>% 
  gather() %>% 
  group_by(key) %>% 

  mode_hdi(value)

# A tibble: 2 x 7
  key    value .lower .upper .width .point .interval
  <chr>  <dbl>  <dbl>  <dbl>  <dbl> <chr>  <chr>    
1 mu    155.   154.   155.     0.95 mode   hdi      
2 sigma   7.74   7.18   8.34   0.95 mode   hdi

Modèle Gaussien sans prédicteur

Taille d'échantillon et normalité

Le prior sur sigma n'est pas une distribution normale mais uniforme sur un intervalle.
Or on vient de voir que la distribution marginale postérieure semble normale

On peut se demander si c'est toujours le cas ???

Considérons le cas où le nombre de mesures est faible… ici 20!

set.seed(4341)
d <- read.csv("Data/Howell1.csv")
d2 <- d %>% filter(age >= 18)
d3 <- sample(d2$height, size = 20)
d3

 [1] 161.925 147.955 156.845 161.290 163.830 161.290 153.035 159.385
 [9] 139.700 154.305 146.050 154.940 149.860 152.400 153.035 160.655
[17] 161.290 148.590 144.780 154.305

Modèle Gaussien sans prédicteur

Taille d'échantillon et normalité

1- définition de la grille

n <- 200
# note we've redefined the ranges of `mu` and `sigma`
d_grid <-
  tibble(
    mu    = seq(from = 150, to = 170, length.out = n),
    sigma = seq(from = 4,   to = 20,  length.out = n)) %>% 
  expand(mu, sigma)

2- définition de la fonction log_likelihood

grid_function <- function(mu, sigma){
  dnorm(d3, mean = mu, sd = sigma, log = T) %>% 
    sum()
}

3- Calcul de la distribution postérieure

d_grid <- d_grid %>% 
  mutate(log_likelihood = map2_dbl(mu, sigma, grid_function)) %>% 

  mutate(prior_mu       = dnorm(mu,    mean = 178, sd  = 20, log = T),
         prior_sigma    = dunif(sigma, min  = 0,   max = 50, log = T)) %>% 

  mutate(product        = log_likelihood + prior_mu + prior_sigma) %>% 

  mutate(probability    = exp(product - max(product)))

4- Échantillonnage sur la distribution postérieure

set.seed(4341)
d_grid_samples <- d_grid %>% 
  sample_n(size = 1e4, replace = T, weight = probability)

Modèle Gaussien sans prédicteur

Taille d'échantillon et normalité

d_grid_samples %>% 
  ggplot(aes(x = mu, y = sigma)) + 
  geom_point(size = .9, alpha = 1/15) +
  scale_fill_viridis_c() +
  coord_cartesian(xlim = range(d_grid_samples$mu),
                  ylim = range(d_grid_samples$sigma)) +
  labs(x = expression(mu[samples]),
       y = expression(sigma[samples])) +
  theme_bw(base_size = 30) +
  theme(panel.grid = element_blank())

Résultat

plot of chunk unnamed-chunk-33

Modèle Gaussien sans prédicteur

Les distributions marginales pour \( ~\mu \) et \( ~\sigma \).

d_grid_samples %>% 
  select(mu, sigma) %>% 
  gather() %>% 

  ggplot(aes(x = value)) + 
   geom_density(colour = "steelblue", size = 1) +
  scale_y_continuous(NULL, breaks = NULL) +
  xlab(NULL) +
  theme_bw(base_size = 30) +
  theme(panel.grid = element_blank()) +
  facet_wrap(~key, scales= "free")

La distribution sur sigma n'est plus vraiment Gaussienne…

Résultat

plot of chunk unnamed-chunk-35

Modèle Gaussien sans prédicteur

Distribution postérieure : méthode par échantillonnage

Le modèle que l'on cherche à fitter sur les données

\[ \begin{align} h_{i} &\sim \mathrm{Normal}(\mu_,\sigma) \\ \mu &\sim \mathrm{Normal}(178, 20) \\ \sigma &\sim \mathrm{Uniform}(0, 50) \end{align} \]

Avec le package brms

Le code correspondant au modèle de gauche.

fit3.1 <- brm(
  data = d2, 
  family = gaussian,
  height ~ 1,
  prior = c(
    prior(normal(178, 20), class = Intercept),
    prior(uniform(0, 50), class = sigma)
    ),
  iter = 31000, warmup = 30000, chains = 4, cores = 4
  )

La fonction qui fit le modèle sur les données est la fonction brm(…)

Le code qui permet d'afficher graphiquement le résultat

plot(fit3.1, theme = theme_bw(base_size = 30))

Modèle Gaussien sans prédicteur

Estimation des paramètres

Le résultat de notre premier modèle en brms

plot(fit3.1, theme = theme_bw(base_size = 30))

plot of chunk unnamed-chunk-39

Modèle Gaussien sans prédicteur

Estimation des paramètres

Le résultat de notre premier modèle en brms.
Description complète

 Family: gaussian 
  Links: mu = identity; sigma = identity 
Formula: height ~ 1 
   Data: d2 (Number of observations: 352) 
Samples: 4 chains, each with iter = 31000; warmup = 30000; thin = 1;
         total post-warmup samples = 4000

Population-Level Effects: 
          Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
Intercept   154.61      0.42   153.80   155.43       4165 1.00

Family Specific Parameters: 
      Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
sigma     7.76      0.29     7.21     8.34        817 1.00

Samples were drawn using sampling(NUTS). For each parameter, Eff.Sample 
is a crude measure of effective sample size, and Rhat is the potential 
scale reduction factor on split chains (at convergence, Rhat = 1).

Description succinte baséée sur la moyenne et la variance (defaut)

posterior_summary(
  fit3.1, 
  probs = c(0.025,0.975), 
  robust = FALSE,) %>% 
  round(digits = 2)

            Estimate Est.Error     Q2.5    Q97.5
b_Intercept   154.61      0.42   153.80   155.43
sigma           7.76      0.29     7.21     8.34
lp__        -1226.87      0.97 -1229.44 -1225.89

Description basée sur la médiane et la déviation absolue

posterior_summary(
  fit3.1, 
  probs = c(0.025,0.975), 
  robust = TRUE,) %>% 
  round(digits = 2)

            Estimate Est.Error     Q2.5    Q97.5
b_Intercept   154.62      0.43   153.80   155.43
sigma           7.74      0.28     7.21     8.34
lp__        -1226.58      0.75 -1229.44 -1225.89

Planning

Cours

Cours 1: Introduction à l'inférence bayésienne

Cours 2: Modèle beta-binomial

Cours 3 : Modèle linéaire à 1 prédicteur continu

Cours 4 : Modèle linéaire multivarié

Cours 5: Markov Chain Monte Carlo

Cours 6: Modèle linéaire généralisé

Cours 7: Comparaison de modèles

Cours 8: Prédicteurs catégoriels et Interactions

Cours 9: Modèles multi-niveaux

Cours 10: Data Hackaton

Plan

Objectifs
Modèle Gaussien sans prédicteur
Modèle Gaussien avec 1 prédicteur
Régression polynomiale

Modèle Gaussien avec 1 prédicteur

Modèle Gaussien sans prédicteur

plot of chunk unnamed-chunk-43

Modèle Gaussien avec 1 prédicteur

plot of chunk unnamed-chunk-45